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Résumé 

Pour modéliser le repliement compact d'un polymère, on étudie 
la statistique des méandres, définis comme configurations d'un circuit 
automobile traversant n fois une rivière. Grâce à une méthode Monte- 
Carlo adaptée à un ordinateur massivement parallèle, les méandres 
sont simulés jusqu'à la taille n = 400. Le comportement asympto- 
tique du dénombrement et de l'enroulement moyen des méandres sont 
présentés. 

Mots-clef : repliement, pliage, méandres, Monte-Carlo, arbre, parallé- 
lisme 



1 Introduction 

Un des centres d'étude actuel de la mécanique statistique est le repliement 
d'objets géométriques qui fluctuent, comme les membranes (à 2 dimensions) 
ou les polymères (à une dimension). Les exemples en physique et en biolo- 
gie sont nombreux. Les concepts développés pour les phénomènes critiques 
(observés lors de certaines transitions de phase) s'appliquent aussi à ce do- 
maine. En particulier, avec "l'universalité", certaines quantités décrivant le 
comportement critique à grande échelle sont indépendantes des détails à pe- 
tite échelle (forme des monomères, interactions, . . . ). Le physicien théoricien 
cherche alors à représenter un objet complexe (par exemple, un polymère bi- 
ologique) par le modèle, le plus simple possible pour pouvoir faire des calculs. 



mais dont on pense qu'il présente les mêmes propriétés physiques à grande 
échelle. 

Dans cet article, le repliement compact d'un polymère est modélisé par 
une bande de timbres-poste complètement repliée [1]. Ceci est équivalent au 
problème des méandres, résumé par cette simple question : de combien de 
manières M„ un circuit automobile peut-il traverser une rivière en n ponts, 
en contournant éventuellement la source ? 

Une variante [2], oii on interdit le contournement de la source, est équi- 
valente à l'énumération des labyrinthes d'un certain type [3], ainsi qu'au 
16ième problème de Hilbert, c'est à dire l'énumération des ovales de courbes 
planaires algébriques [4]. 

Ce problème a par ailleurs des connections avec la théorie des nœuds, les 
"modèles de matrice" et la "gravité quantique", et même la QCD (théorie 
des interactions nucléaires fortes en physique des particules). 

Les méandres sont donc un modèle de repliement compact d'un objet à 
une dimension, avec des maillons identiques, dont on ne retient que l'ordre 
dans le pliage. Malgré cette simplification extrême, ce problème résiste depuis 
au moins un siècle : par exemple, on ne connaît toujours pas de formule pour 
le nombre de pliages possibles M„, ni même son comportement asymptotique 
pour n grand. 

L'approche présentée ici est purement numérique. Elle a été rendue pos- 
sible en exploitant la puissance des ordinateurs parallèles. Cependant, la 
méthode de "brute force" ne permet pas d'aller bien loin et un algorithme 
Monte-Carlo efficace a été élaboré. La nature même des méandres rend diffi- 
cile la vectorisation des boucles les plus internes du programme. Par contre, 
celui a été parallélisé avec des "gros grains" , c'est à dire en divisant le travail 
en sous-tâches au plus haut niveau possible. 

La Section 2 est consacrée aux définitions. La Section 3 présente un al- 
gorithme de dénombrement exact, lui- aussi parallélisé. La Section 4 décrit la 
méthode Monte-Carlo. Finalement, quelques résultats sont exposés Section 5. 

2 Définitions 

Un méandre de taille n est défini comme une configuration planaire d'une 
boucle qui ne se croise pas elle-même (la route) , croisant une demi-droite (la 
rivière avec sa source) en n points (les ponts). Deux méandres de même taille 
sont considérés comme équivalents si l'un est obtenu à partir de l'autre en 
déformant continûment la route tout en gardant fixés les ponts. C'est une 
équivalence topologique. 





FiG. 1: Les M4 = 4 méandres de taille 4. La route est la boucle auto- évitante. 
La rivière semi-infinie est la demi-droite en trait plein, commençant à la 
source (point noir). La taille est le nombre de ponts. Le nombre d'enroule- 
ments w est le nombre d'arches croisant la demi- droite en pointillé à droite 
de la source. Les deux méandres du haut n'ont pas d'enroulement (w = 0), 
mais les deux du bas ont w = 2. 

On appelle arche chaque section de route entre deux ponts. Un méandre 
de taille n a donc n ponts et n arches. 

Le nombre de méandres différents de taille n est noté M„. Par exemple, 
Ml = 1, M2 = 1, M3 = 2 et M4 = 4. Sur la Fig. 1, les 4 méandres de taille 4 
sont dessinés. 

Le problème des méandres est absolument équivalent à celui du repliement 
d'une bande de timbres-poste. Considérons une bande de n timbres attachée 
à un support. Il y a alors M„_|_i manières de plier cette bande complètement, 
c'est à dire sous forme de pile avec un seul timbre en largeur. En effet, comme 
expliqué Fig. 2, chaque pliage de n timbres correspond à un méandre de taille 
n + 1 et réciproquement. On préfère cependant utiliser dans cet article la 
représentation sous forme de méandre car elle est plus naturelle pour décrire 
la récurrence avec laquelle nous les simulerons. 

Les méandres ont des similitudes avec d'autres problèmes de physique 
statistique, par exemple les marches aléatoires auto-évitantes fermées : un 
méandre est obtenu à partir d'une telle marche en lui superposant une demi- 
droite et en ne gardant que l'aspect topologique, c'est à dire la nature des 
intersections. Aussi, par analogie, on attend que 

n -- C—, 1 

oii les estimations [5, 6] sont R ~ 3.5 et 7 ~ 2. 

Le nombre R, équivalent de la connectivité, peut être interprété comme 
le nombre moyen de possibilités pour ajouter un pont près de la source en ne 
déformant qu'une seule arche, pour un grand méandre. Ainsi, ln{R) apparaît 
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(a) (b) 

FiG. 2: On transforme un pliage de n timbres (ici n = 4) en un méandre de 
taille n + 1 de la manière suivante : (a) dessiner une boucle ovale traversant 
la pile de n timbres et se refermant sur la gauche en traversant le support. Il 
y a donc n + 1 intersections, (b) Mettre à plat la bande de timbres en tirant 
vers la droite son extrémité. On obtient alors un méandre avec n+1 ponts, la 
boucle devenant la route, et la bande de timbres la rivière. La transformation 
inverse donne, à partir de tout méandre, une bande de timbre totalement 
repliée, en déformant la route pour en faire une boucle ovale. 

comme l'entropie par pont ; c'est une quantité intensive, c'est à dire définie 
par unité de volume. Par contre, l'exposant 7 sera sensible aux conditions au 
bord. Par exemple, si l'on change la définition du problème en prenant une 
rivière infinie aux deux bouts, ou bien une rivière avec un affluent, alors la 
connectivité R ne sera pas changée (elle caractérise l'effet du volume), mais 
par contre 7 sera à chaque fois différent. Numériquement, 7 est plus difficile 
à mesurer car il décrit la correction sous-dominante dans l'équation (1). 

Pour un méandre donné, le nombre d'enroulements de la route autour de 
la source est défini comme le nombre minimal d'intersections entre la route et 
une demi-droite (semi-infinie) partant de la source en prolongeant la rivière 
du coté opposé. On peut le voir comme la distance topologique bout-à-bout 
entre la source (extrémité droite de la rivière) et l'infini (extrémité gauche), 
la distance entre deux points étant définie comme le nombre minimal d'arches 
à traverser pour aller d'un point à l'autre, sans franchir la rivière. Pour un 
exemple, voir Fig. 1. 

On définit w„ comme la moyenne du nombre d'enroulements sur tous les 
Mn méandres de taille n. Par analogie avec les marches aléatoires, on s'attend 
a 

n— >oo i, /„x 

Wn ^ n , (2) 

avec < z/ < 1. 

On peut étudier une variante des méandres oii le nombre de routes est 
laissé libre ; les routes ne doivent pas se croiser, ni elles-mêmes, ni entre elles 
(mais peuvent être emboîtées) et le nombre total de ponts sur la rivière est 



toujours fixé a n. Ce problème est alors exactement équivalent [7] à celui d'une 
marche aléatoire de 2n pas sur une demi-droite, le marcheur commençant et 
finissant à l'extrémité. Chaque marche correspond alors à une configuration 
de méandres et réciproquement. Il y a alors C„ marches possibles de ce type, 
où les Cn sont les nombres de Catalan (2n)!/(n + l)!/n!. On en déduit, pour 
ce modèle soluble, que i? = 4, 7 = 3/2 et z/ = 1/2 qui est l'exposant du 
mouvement Brownien. Mais en imposant qu'il n'y ait qu'une seule route, la 
nature du problème change complètement et les résultats exacts sont rares. 
En particulier, les valeurs ci-dessus changent et n'ont pas encore été calculées 
exactement. 

On peut aussi définir d'autres quantités comme la forme moyenne des 
méandres ou leur aires. Plus de détails sont donnés dans la Réf. [8]. 

3 Dénombrement exact 

Pour construire les méandres de taille n-|-l à partir des méandres de taille 
n de façon systématique, plusieurs manières sont possibles. Malheureusement, 
aucune d'entre elles ne donne de récurrence directement entre le nombre 
total de méandres M„_|_i, et M„ : il faut étudier séparément chacun des M„ 
méandres. Dans cette section, nous décrivons une méthode [6, 7], qui a été 
programmée sur un ordinateur parallèle, le Cray-T3D. Comme elle est a la 
base de l'algorithme Monte-Carlo, nous la présentons en détail. 

Cette méthode consiste à ajouter un pont sur la partie la plus à gauche 
de la rivière et à modifier la route pour la faire passer par ce nouveau pont. 
Pour que ce changement soit minimal, seule une arche extérieure est modifiée 
(une arche est dite extérieure quand aucune autre ne l'entoure). 

Prenons un méandre de taille n (le père) et choisissons (ou marquons) une 
de ses arches extérieures. Par le processus décrit Fig. 3, un méandre de taille 
n + 1 (le fils) est construit. Le père a autant de fils que d'arches extérieures. 
En inversant le processus, chaque méandre de taille n + 1 a un et un seul 
père. C'est donc une correspondance un-à-un entre, d'une part les méandres 
de taille n + 1, et d'autre part, ceux de taille n avec une arche extérieure 
marquée. 

Le point de départ de la récurrence est l'unique méandre de taille L En 
appliquant n — 1 fois le processus, tous les méandres de taille n peuvent 
être construits. Comme décrit Fig. 4, l'ensemble des méandres s'organise 
comme un arbre. La racine, au niveau 1, est le méandre de départ n = 1. 
Chaque branche entre un nœud au niveau n et un autre nœud au niveau 
n + 1 représente une relation entre un père de taille n et un de ses fils de 
taille n + 1. A l'exception de la racine n = 1, chaque méandre (ou nœud) a 




FiG. 3: Un méandre (fils) de taille n + 1 est construit à partir d'un père 
de taille n avec une arche extérieure marquée de la manière suivante : (a) 
Ajouter un pont sur la partie gauche de la rivière. Couper l'arche extérieure 
marquée et tirer sur ses deux extrémités libres, (h) Refermer cette arche du 
coté opposé en passant par le nouveau pont (et éventuellement en contournant 
la source par la droite). Le processus est inversible : (b) Ouvrir la route à 
l'emplacement du pont le plus à gauche (a) Prendre les deux extrémités libres 
et les refermer du coté opposé pour former une arche extérieure. 

plusieurs arches extérieures, donc plusieurs fils (ou branches). Leur nombre 
dépend de la forme précise du méandre père et varie entre 2 et n/2 + 1. Il 
s'agit donc d'un arbre déterministe mais dont le nombre de branches dépend 
du nœud considéré. 

Si l'on veut énumérer exactement les M„ méandres de taille n, la seule 
méthode que nous connaissons est de construire cet arbre jusqu'au niveau n. 
En particulier, nous n'avons pas trouvé de relation de récurrence directement 
entre les M„. Le nombre de fils pour chaque méandre a une distribution qui 
semble erratique et le seul moyen de la connaître est l'examen de ses arches. 
Aussi, pour calculer Mn, le travail à fournir est proportionnel aux M„, qui 
croissent exponentiellement : les limites de l'ordinateur sont vite atteintes. 

Pour cela, on coupe donc l'arbre au niveau n. On appelle alors feuilles les 
nœuds de niveau n (qui représentent des méandres de taille n). L'algorithme, 
schématisé Fig. 5, consiste à visiter toutes les feuilles de la gauche vers la 
droite, à la manière d'un "écureuil méticuleux" . Les règles sont les suivantes. 
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FiG. 4: L'arbre des méandres jusqu'à n = 5. Pour n > A, seule une moitié 
est dessinée, l'autre étant obtenue par la symétrie haut-bas. Ainsi, Mi = 1, 
M2 = l, M3 = 2, M4 = 4 et M5 = 10. Chaque flèche représente un processus 
"père-fils". 

(a) L'écureuil démarre à la racine, (b) Lorsque l'écureuil se trouve sur un 
nœud intermédiaire (qui n'est pas une feuille), il grimpe dans la branche 
la plus à gauche qu'il n'a pas encore visitée. S'il a déjà visité toutes les 
branches, il descend d'un niveau, (c) Lorsque l'écureuil se trouve sur une 
feuille, il descend d'un niveau. 

Le lecteur peut se convaincre que ces règles décrivent bien une visite 
complète de l'arbre. Bien sur, lorsque l'écureuil se trouve sur un nœud, il 
mesure des quantités intéressantes, par exemple le nombre d'enroulements. 
Ces mesures sont cumulées et traitées à la fin de l'énumération. 

D'un point de vue Fortran, plusieurs représentations des méandres sont 
possibles. Dans celle que nous avons finalement utilisée, chaque arche est 
codée par les numéros des deux ponts qu'elle relie. Il faut cependant dis- 
tinguer les deux rives. Aussi, les ponts, vus de la rive supérieure, sont numérotés| 
négativement de — n + 1 à 0, de la gauche vers la droite. Vus de la rive 
inférieure, ils sont numérotés positivement de 1 à n, de la droite vers la 
gauche. Le pont le plus proche de la source porte les numéros et 1 ; le plus 
éloigné les numéros —n + 1 et n. Ainsi, le système des n arches est codé par 
le tableau d'entiers A{—n + 1 : n) 011 A{i) est le numéro du pont relié par 
une arche au pont i. Voir Fig. 6 pour un exemple. 

Ce codage est loin d'être optimal pour l'utilisation de la mémoire. On 
pourrait par exemple se contenter d'un tableau de 2n bits, 011 le bit i vaut 1 




FiG. 5: Algorithme de l'écureuil méticuleux : il visite toutes les feuilles de la 
gauche vers la droite. 




FiG. 6: 

-V- 



Pour ce méandre de taille 4, l'arche (—1,4) est représentée par 
= A et A{A) = -1. Au total, A{-3 : 4) = (-2, -3, A, 3, 2, 1, 0, 



(resp. 0) si le pont i est relié à un pont de numéro supérieur (resp. inférieur). 
Cependant ce codage permet de casser et de fusionner des arches efficace- 
ment. En effet, lorsque l'écureuil monte, une arche {j,A{j)) est cassée et 
donne deux arches limitées par le nouveau pont : {—n,j) et {A{j),n + 1). 
Lorsque l'écureuil redescend, les deux arches {—n + l,A{—n + l)) et {A{n),n) 
fusionnent pour donner l'arche {A{—n+ l),A{n)), l'arche suivante à casser 
étant éventuellement celle commençant k j = A{n) + 1. 

Le programme Fortran, Table 1, compte donc les méandres jusqu'à n = 
nmax donné. Il est possible d'utiliser la symétrie haut-bas pour diviser le 
travail par deux. On peut aussi se contenter de ne construire les méandres 
que jusqu'à la taille n — 1 et le décompte de leurs arches externes donne M„. 
On peut gagner k étages supplémentaires, en utilisant des formules, expo- 
nentiellement compliquées avec k, faisant intervenir pour chaque méandre de 
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Integer, 


Parameter : : nmax = 20 


Integer : 


: A(-nmax+l:nmax) = 


Integer : 


: n = 1 


Integer : 


: j 


Integer : 


: M (nmax) = 


A(0) = 1 


Arbre : Do 


M(n) = M(n) + 1 


j = -n + 


1 



Do While (n==nmax .or. j==n+l) 

A(A(-n+l)) = A(n) 

A(A(n)) = A(-n+l) 

j = A(n) + 1 

n = n - 1 

If (n == 0) Exit Arbre 
Enddo 



hauteur maximale 
représentation en arches 
hauteur de l'écureuil 
branche a visiter 
compte les méandres de 
taille n 
méandre de départ 

parcours de l'arbre 
visite d'un nouveau noeud 
arche la plus a gauche 

Descente : 
on fusionne les deux 

arches externes 
arche suivante a casser 

fin du parcours 



A(A(j)) = n+1 
A(n+1) = A(j) 
A(j) = -n 
A(-n) = j 
n = n + 1 
Enddo Arbre 



! Montée : 

! l'arche (j) est cassée 



Print '(i3,il5)', (n, M(n) , n = 1, nmax) 
End 



Tab. 1: Programme Fortran qui compte les méandres jusqu'à n 



nmax. 



PO 



PI 



P2 



PO 



PI 




FiG. 7: Parallélisation de l'algorithme de l'écureuil méticuleux : dans une 
phase initiale très courte, tous les processeurs construisent l'arbre jusqu'à une 
hauteur Ui. Les méandres obtenus sont alors pris comme nouveaux points de 
départ et le travail est réparti de manière cyclique entre tous les processeurs, 
notés ICI PO, PI et P2. 

taille n — k, les nombres de ses arches externes, de ses arches de profondeur 
2, . . . , jusqu'à la profondeur k. Le choix optimal est fc = 4 ou 5. 

Ce programme est clairement anti-vectoriel. Pour donner un ordre de 
grandeur, avec n = 20, M20 = 102511418 est obtenu en 31 secondes sur 
un Sun avec UltraSparc à 360 MHz. Il met 70 secondes sur une machine 
vectorielle : le Fujitsu VPP-300 du Cea-Grenoble. Pour étudier les n les 
plus grands possible avec cet algorithme (qui est le seul dont on dispose), il 
ne reste plus que la voie du parallélisme massif, 011 l'ordinateur dispose de 
nombreux processeurs calculant indépendamment et interconnectés par un 
réseau rapide. 

Pour cela, une taille intermédiaire ni est choisie. Dans une première phase, 
très courte, l'arbre est construit, sur chaque processeur, jusqu'au niveau ni. 
On obtient donc M„^ "petits" méandres, par exemple Mu = 4210 avec ni = 
11. 

Dans une seconde phase, longue et massivement parallèle (voir Fig. 7), 
chacun de ces petits méandres est considéré comme la racine d'un arbre, 
qui est alors un sous-arbre de l'arbre global. Chaque sous-arbre est traité 
indépendamment des autres. Sur une machine parallèle, le travail est réparti 
entre tous les processeurs, chacun traitant, à peu prés, le même nombre de 
sous-arbres. A la fin, les mesures sont collectées, sommées et traitées par un 
seul processeur. 

Pour un sous-arbre donné, le travail à faire est proportionnel à son nom- 
bre de nœuds : il fluctue sensiblement d'un sous-arbre à l'autre. Mais comme 
chaque processeur a plusieurs dizaines de sous-arbres à traiter, on constate 
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n 


M„ 


n 


M„ 


1 


1 


16 


1053874 


2 


1 


17 


3328188 


3 


2 


18 


10274466 


4 


4 


19 


32786630 


5 


10 


20 


102511418 


6 


24 


21 


329903058 


7 


66 


22 


1042277722 


8 


174 


23 


3377919260 


9 


504 


24 


10765024432 


10 


1406 


25 


35095839848 


11 


4210 


26 


112670468128 


12 


12198 


27 


369192702554 


13 


37378 


28 


1192724674590 


14 


111278 


29 


3925446804750 


15 


346846 







Tab. 2: Les nombres Mn de méandres de taille n pour n < 29, obtenus par 
énumération exacte sur ordinateur parallèle. 

que les fluctuations se moyennent et l'équilibrage des charges entre pro- 
cesseurs est réalisé à quelques % près. Il n'est donc pas nécessaire d'utiliser 
de méthode plus sophistiquée. 

Pour programmer cela sur un Cray-T3D, les directives de compilation 
Craft, fournies par le constructeur, ont été utilisées, car l'algorithme s'y 
prêtait bien. Cette extension de Fortran est du type "data-parallel" et a 
été remplacée depuis par Hpf. Une directive barrier synchronise tous les pro- 
cesseurs avant la collecte des résultats, qui se fait par le biais d'un tableau 
partagé, oii chaque processeur y écrit le résultat de son dénombrement. Puis, 
grâce à des directives master et end master, seul un processeur fait et traite 
les sommes. La difficulté principale n'a donc pas été le codage en soi de 
l'algorithme dans le langage parallèle, mais bien la mise au point d'une par- 
allélisation efficace. 

Les résultat jusqu'à n = 29, donnés Table 2, ont été obtenus en 163 heures 
réparties sur 128 processeurs du Cray-T3D du Cea-Grenoble en 1995. On voit 
que les M„ croissent exponentiellement. Après extrapolation pour n = oo, 
on en déduit que R = 3.50(1) et 7 ~ 2. 

Si la puissance des ordinateurs continue de croître exponentiellement au 
fil des années, comme le temps de calcul est proportionnel aux M„ qui crois- 
sent eux-aussi exponentiellement, le mieux qu'on puisse espérer avec cette 
méthode de dénombrement exact est un gain d'une nouvelle taille (de n à 
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n + 1) tous les deux ans environ. 

4 Méthode Monte-Carlo 

Nous avons vu dans la section précédente que la méthode d'énumération 
exacte était vite limitée car le nombre M„ de méandres à construire croît 
exponentiellement avec la taille n. L'idée générale des méthodes Monte-Carlo 
est de remplacer l'analyse complète de tous les cas possibles par la sélection 
aléatoire d'un sous-ensemble. Celui sera représentatif si la moyenne obtenue 
avec ces échantillons est proche du vrai résultat (qui est inconnu). Cette 
moyenne est aléatoire et fluctue d'une expérience numérique à l'autre, ces 
fluctuations étant alors empiriquement mesurables. 

Pour le problème des méandres, comme on ne sait construire efficacement 
un méandre de taille n + 1 qu'à partir d'un méandre de taille n, nous allons 
nous servir de l'arbre décrit plus haut. 

Nous allons d'abord introduire la méthode qui utilise un écureuil Monte- 
Carlo. Malheureusement cette méthode n'est pas efficace car ses fluctuations 
statistiques croissent très vite avec n. La parade consiste alors à utiliser une 
grande population d'écureuils, méthode qui, de plus, est parallélisable. 

4.1 Un écureuil Monte-Carlo 

La vie de l'écureuil Monte-Carlo commence à la racine de l'arbre (à n = 
1). Il grimpe de manière aléatoire dans l'arbre. A chaque niveau n, il se 
trouve sur un nœud de l'arbre (qui représente un méandre de taille n) et 
fait un certain nombre de mesures. Il grimpe au niveau n + 1 en choisissant 
au hasard de manière uniforme l'une des 6„ branches partant de ce nœud. 
L'écureuil s'arrête à une certaine hauteur Umax fixée à l'avance. 

Comme représentée Fig. 8, la probabilité que l'écureuil atteigne au niveau 
n un nœud donné n'est pas uniforme. Elle vaut l/g„ où 

qn = h-h2---K~i (3) 

est le produit des nombres de branches hi que l'écureuil a rencontrés à chaque 
nœud lors de son ascension, et qu'il peut donc calculer. 

Ce processus, c'est à dire une ascension Monte-Carlo d'un seul écureuil, 
est considéré comme une simulation. De nombreuses simulations indépen- 
dantes sont réalisées et les mesures sont moyennées. Pour cela, il faut les 
pondérer avec le poids g„. Ainsi, chaque méandre contribue au résultat avec 
probabilité l/g„ et, le cas échéant, le poids g^, donc de manière uniforme. 
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FiG. 8: La probabilité (nombres du haut) d'atteindre une feuille donnée est 
le produit des probabilités d'embranchement (nombres intermédiaires) ren- 
contrées à chaque noeud. Elle n'est pas uniforme. 

En particulier, l'espérance mathématique (g„), définie comme la moyenne sur 
toutes les simulations possibles, vaut M„. 

Malheureusement, cette méthode ne marche pas en pratique pour n grand. 
En effet, bien que la loi de probabilité de chaque bi soit régulière, la loi des g„, 
produit d'un grand nombre de bi, n'est pas auto-moyennante : dans la limite 
n grand, avec probabilité 1, les g„ obtenus par simulation sont bien plus petits 
que l'espérance mathématique (g„). En effet, le théorème de la limite centrale 
s'applique à lng„, qui est auto-moyennant car somme des ln6j, mais pas à g„- 
Plus précisément, des événements exponentiellement rares avec n contribuent 
de manière exponentiellement grande à (ç„). Dans une simulation, la moyenne 
des Qn observés sera dominée par ces événements rares et les fluctuations 
seront larges. 11 faut donc un nombre exponentiellement grand (avec n) de 
simulations. Dans la pratique, on ne peut pas dépasser n ~ 30 ou 35. 

4.2 Méthode Monte-Carlo multi-écureuil 

La méthode précédente est généralisée en utilisant une population de S 
écureuils, oii S est un paramètre fixé à l'avance. Il est plus simple de le choisir 
parmi la liste des M„ : S = Mng. Ainsi, au départ, chaque nœud au niveau hq 
est occupé par un écureuil. Dans ce travail, Uq = 17 et la symétrie haut-bas 
a été utilisée pour réduire la population h S = Mn/2 = 1 664 094 écureuils. 

Les S écureuils évoluent, de manière synchrone, du niveau n au niveau n + 
1 par le processus suivant, schématisé Fig. 9. Considérons l'écureuil d'index i 
sur un nœud au niveau n, connecté à bi nœuds du niveau n+1. Il se reproduit 
et donne bi fils, chacun d'entre eux vivant sur un de ces bi nœuds supérieurs. 
Tous les écureuils se reproduisent en même temps et le nombre total de fils 
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FiG. 9: La méthode Monte-Carlo multi- écureuil est ici schématisée avec S = 
2 écureuils (El et E2). Leurs S' = 5 fils (FI à F5) sont énumérés. Pour 
garder S constant, ils sont décimés et seulement S = 2 fils, choisis au hasard, 
survivent. Ce processus est itéré jusqu'à atteindre la hauteur voulue l'arbre. 

S' = Y.i=i bi est calculé. Le taux de reproduction moyen _B„ = S' /S est une 
estimation du rapport M„_|_i/M„. Pour empêcher la croissance exponentielle 
de la population (et donc des besoins en mémoire et en Cpu), celle-ci est 
gardée constante : seulement S parmi les S' fils survivent. Le choix est fait au 
hasard de manière uniforme. La probabilité de survie d'un fils est donc 1/Bn- 
Cette décimation est l'étape Monte-Carlo de l'algorithme. Puis, le processus 
s'applique à nouveau aux fils ; il est itéré jusqu'à atteindre une hauteur Umax 
fixée à l'avance. Cela constitue une simulation. De nombreuses simulations 
indépendantes sont réalisées : les moyennes et leurs barres d'erreur associées 
sont alors calculées. 

On voit que le cas particulier S = 1 n'est rien d'autre que la méthode 
précédente à un écureuil. A l'opposé, la limite S* = oo est équivalente à la 
méthode d'énumération exacte car on ne supprime pas d'écureuil. 

Mais quelque soit S, comme pour S* = 1, la probabilité qu'un nœud soit 
visité par un écureuil n'est pas uniforme et il y a toujours un biais. Par 
exemple, les nœuds qui ont peu de frères, de cousins, . . . , ont plus de chance 
d'être visités. On peut montrer [8] que ce biais est exactement corrigé si les 
mesures sont pondérées, pour chaque simulation, par 

Qn = Bna . . . Bn-2-Bn-lj (4) 

produit des facteurs de décimation. 

A priori, cette méthode souffre des mêmes défauts que la précédente, car le 
poids qn est le produit de nombreux Bi aléatoires, donc non auto-moyennant 
quand n devient grand. Mais, l'amélioration capitale est que la distribution 
de Bi devient très piquée autour de Mj+i/Mj quand la population S est 
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grande. En effet, Bi étant le nombre de fils en moyenne pour S écureuils, ses 
fluctuations sont faibles et d'ordre 0(l/viS). Pour le dire autrement, dans 
l'expression de g„, la moyenne sur les S écureuils est faite au niveau des 
Bi, donc avant leur produit. Les fluctuations de g„ croissent certes toujours 
exponentiellement, mais avec un taux très faible si S est grand. Par exemple, 
pour S* = 1 664 094 et n = 400, la largeur a de la distribution des g„ est de 
12% avec quelques événements rares allant jusqu'à trois fois la moyenne. 

Comment choisir S et n7 Naïvement, on a envie de simuler les méandres 
de taille n la plus grande possible. Mais, pour réaliser A^^ simulations indépen- 
dantes avec S écureuils jusqu'à la taille n, le besoin de mémoire d'ordina- 
teur est d'ordre 0{n.S) et celui de temps Cpu d'ordre O(n^.S'.A^s). Avec S 
suffisamment grand, les fluctuations statistiques sont presque gaussiennes et 
d'ordre 0{l/\/S.Ns). Comme n est toujours limité, nous extrapolerons pour 
étudier la limite n infini. L'allocation de temps d'ordinateur étant fixée par 
ailleurs, nous nous sommes limités à n < 400, pour bien extrapoler avec 
une meilleure statistique. Finalement, à produit S.Ng constant, il vaut mieux 
choisir S le plus grand possible permis par la mémoire de l'ordinateur, pour 
diminuer le problème des fluctuations rares et grandes de g„. 

L'algorithme a été ici décrit en terme d "écureuils" . Mais ceux-ci repré- 
sentent des méandres par l'intermédiaire d'un tableau d'entiers A{—n+l : n). 
En particulier, l'opération qui transforme un écureuil en son fils, consiste à 
casser une arche externe, donc à manipuler ce tableau. Cette opération n'est 
pas vectorisable et s'avère être la partie qui consomme le plus de Cpu dans 
l'algorithme. Comme dans toute simulation Monte-Carlo, la qualité statis- 
tique des résultats est limitée par le temps Cpu. Il est alors capital de par- 
alléliser l'algorithme. 

La parallélisation la plus simple consiste à réaliser des simulations totale- 
ment indépendantes sur chacun des processeurs. Il suffit d'initialiser diffé- 
remment le générateur de nombres aléatoires d'un processeur à l'autre, et 
d'organiser la collecte des mesures effectuées lors des simulations pour les 
traiter ultérieurement. Ce type de parallélisation n'utilise que très peu le 
réseau car rien n'est échangé entre les processeurs. Cela peut même se con- 
cevoir sur un réseau hétérogène de stations de travail, de manière asynchrone. 

Mais on a vu ci-dessus qu'il était très important de simuler la plus grande 
population possible, plutôt que de multiplier les simulations avec une popula- 
tion plus petite, la limite étant donnée par la mémoire disponible. Pour cela, 
tous les processeurs doivent travailler de concert sur la même population. 
L'algorithme utilisé est décrit ci-après. 

Dans une courte phase d'initialisation, les S écureuils sont construits au 
niveau Uq et sont répartis de manière équilibrée parmi les processeurs. Chaque 
processeur ne s'occupe que de son groupe d'écureuils. 
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Lors d'une itération principale, oii la population passe du niveau n au 
niveau n + 1, chaque processeur énumère les fils de ses écureuils, puis réalise 
et somme ses mesures. Les sommes globales sont ensuite faites sur l'ensemble 
des processeurs et écrites dans un fichier par un seul processeur pour être 
analysées ultérieurement. La facteur de décimation i?„ doit être déterminé de 
manière globale, et non pas localement sur chaque processeur. Un processeur, 
déclaré maître, analyse le nombre de fils Fp obtenus par chaque processeur 
p, calcule Bn et fixe en retour le nombre Fp/ Bn de fils que chaque processeur 
doit garder. La population globale est donc maintenue exactement égale à S. 
Cette phase demande une synchronisation entre processeurs, peu de calculs 
et l'échange de petits messages. Puis vient la décimation proprement dite, 
réalisée de manière indépendante sur chacun des processeurs. C'est l'étape qui 
consomme le plus de Cpu car les fils survivants sont effectivement construits 
à ce moment- là. 

Il est alors nécessaire d'équilibrer les populations entre processeurs. En 
effet, le nombre de fils fiuctue. Lors de la synchronisation, tous les processeurs 
doivent attendre celui qui a le plus d'écureuils, car il effectue plus de travail. 
Sans y remédier, ce problème empirerait au fil des itérations. Aussi, à la 
fin de chaque itération, le processeur maître recense les processeurs dont 
la population dépasse de 5 % la moyenne et leur fait transférer la partie 
excédentaire aux processeurs les moins chargés. Ainsi, les temps d'attente 
sur la barrière de synchronisation sont en moyenne de 5 %. En contrepartie, 
le réseau transfère 5 % de la mémoire de quelques processeurs. 

Ce processus constitue une simulation. Plusieurs simulations sont ef- 
fectuées de manières indépendantes, en un ou plusieurs "runs". Les fichiers 
de mesures sont ensuite traités par un autre programme. 

Cette parallélisation a été programmée sur le Cray-T3E du Cea-Grenoble, 
équipé de processeurs Dec-alpha à 375 MHz. Il faut 13 Goctets de mémoire 
pour une population de S* = 1 664 094 écureuils, soit 128 processeurs de 
128 Moctets. Pour réaliser 8192 simulations, 8 jours ont été nécessaires, soit 
24000 heures de Cpu. Pour paralléliser, la bibliothèque Shmem, fournie par 
le constructeur, a été utilisée. Elle est du type "passage de message" , comme 
PVM ou MPI. Son intérêt principal est sa simplicité et son efficacité. Par 
contre, elle n'est pas portable. 

Le rendement de la parallélisation est proche des 95 % théoriques. Aucune 
dégradation significative n'a été vue, même avec 128 processeurs. En effet, 
les échanges entre processeurs sont rares et de volume moyen : ils ne posent 
aucun problème pour le réseau du T3E, qui est un point fort de cette machine. 
La performance du programme est plutôt sensible aux problèmes d'accès au 
cache et à la mémoire locale ; un gain important a été obtenu, lors des copies 
de tableaux dans une même mémoire locale, grâce aux registres "E", en 
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court-circuitant le Cpu et son cache. 

On voit que le travail scalaire par processeur est proportionnel à son 
nombre d'écureuils ; par contre la taille des messages est proportionnelle aux 
fluctuations, donc à la racine carrée. En théorie, la "scalabilité" est par- 
faite, à condition d'augmenter la population globale — ce qui est justement 
intéressant — quand on augmente le nombre de processeurs, pour que la 
population par processeur ne s'appauvrisse pas. 

5 Quelques résultats 

On s'intéresse tout d'abord au comportement asymptotique du nombre 
de méandres M„ ~ cR"" /rC pour n grand. L'entropie In R peut être estimée 
avec 

i»4lnf^V (5) 



OÙ n saute de 2 en 2 pour diminuer des effets sensibles de parité entre les n 
pairs et impairs. Comme on s'attend à L„ ~ Ini? — 7/n pour n grand, en 
traçant y = Ln en fonction de a; = 1/n, on pourra estimer Ini? (limite quand 
X tend vers 0) et 7 (pente asymptotique). 

Sur la Fig. 10, l'estimation Monte-Carlo de L„ — ln3.5 + 2/n est tracée en 
fonction de 1/n pour n allant de 50 à 400. La fonction linéaire y = In 3.5 — 2x 
a été arbitrairement soustraite pour réduire l'amplitude de variation de y ; 
les quantités intéressantes 2 — 7 (pente résiduelle), ln(i?/3.5) (limite quand 
X tend vers 0) et la courbure (déviation au comportement linéaire espéré) 
sont ainsi rendues plus visibles. Malgré cette amplification, la courbure reste 
faible et rend raisonnable une extrapolation linéaire, qui donne 

iî = 3.5019(2) et 7 = 2.056(10), (6) 

en contradiction avec les conjectures 7/2 et 2. 

U exposant d'enroulement u, défini par le comportement asymptotique du 
nombre moyen d'enroulements Wn ~ n^, est présenté Fig. 11. En traçant 
ln{wn) en fonction de Inn, la pente asymptotique est une mesure de z/. Nous 
préférons utiliser ln(ti;„ + 1) au lieu de \n{wn), car on remarque que w„ + 1 est 
moins sensible que w„ aux effets de taille finie. Comme la question principale 
est de savoir si z/ = 1/2 ou non, la fonction linéaire y = x/2 a été soustraite. 
Ainsi, la variation de y est réduite ; z/ — 1/2 et la courbure sont plus visibles. 
On voit que cette dernière reste petite et une extrapolation linéaire donne 

z/ = 0.518(2), (7) 
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FiG. 10: Estimation Monte-Carlo de Ln — In 3.5 -|- 2/n en fonction de 1/n 
pour n allant de 50 à 400. La limite quand x tend vers est ln(i?/3.5), et la 
pente (négative) est2~'~f. Les barres d'erreur ne sont pas dessinées car elles 
sont toujours inférieures à 10~^, donc plus petites que les symboles. Un effet 
de parité entre les n pairs et impairs est visible. Une extrapolation linéaire 
donne R = 3.5019(2) et 7 = 2.056(10). 
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FiG. 11: Estimation Monte-Carlo de In(-U7„-|- 1) — 1/2 Inn en fonction de Inn, 
pour n entre 8 et 400. La pente est de 0.018 et mesure v — 1/2. Les barres 
d'erreur ne sont pas dessinées car inférieures à 3.10^^. 
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incompatible avec la valeur brownienne 1/2. 

D'autres résultats sont présentés dans la Réf [8], concernant la loi de 
probabilité du nombre d'enroulements et la forme moyenne des méandres, en 
particulier près de la source et de l'extrémité opposée. 

6 Conclusion 

Dans ce travail, nous avons simulé le repliement compact d'un objet uni- 
dimensionnel, une bande de n timbres-poste, modélisée par les méandres. 
Jusqu'à présent, on ne connaît pas de formule donnant le nombre de pliages 
possibles M„ en fonction de n, et l'ordinateur est donc une méthode d'inves- 
tigation privilégiée. Les algorithmes, que nous avons développés et utilisés, 
reposent sur une représentation de l'ensemble de tous les pliages possibles en 
un arbre de forme irrégulière. 

Avec un algorithme de visite systématique de tous les nœuds de l'arbre 
jusqu'à une hauteur fixée n, on peut mesurer exactement toutes les pro- 
priétés. Mais la nature même de l'arbre empêche la vectorisation, car il n'y 
a pas de longue boucle régulière. Par contre, la possibilité de le découper en 
sous-arbres rend la parallélisation efficace sur une machine à grand nombre 
de processeurs. Mais le travail croît exponentiellement avec n et même un 
ordinateur puissant ne permet pas de dépasser n ~ 30. 

Avec un algorithme Monte-Carlo, on peut étudier des méandres plus 
grands. Mais les fluctuations statistiques, inhérentes à ce type de méthode, 
sont ici particulièrement défavorables. Avec la mise au point de la méthode 
multi-écureuil, les moyennes sont réalisées en deux étapes : la première, entre 
écureuils, a des fluctuations gaussiennes, le cas le plus favorable. On cherche 
donc à simuler simultanément le plus grand nombre possible d'écureuils, ce 
qui est principalement limité par la mémoire de l'ordinateur. En répartissant 
la population sur de nombreux processeurs, la parallélisation permet alors de 
bénéficier de beaucoup plus de mémoire. On obtient aussi une forte puissance 
en Cpu scalaire, nécessaire pour que les résultats significatifs ne soient pas 
noyés dans les fiuctuations statistiques, ce qui a permis alors des simulations 
jusqu'à la taille n = 400. 

On voit donc que les études décrites dans cet article doivent beaucoup 
aux ordinateurs parallèles, les ordinateurs vectoriels étant inefficaces pour 
ces algorithmes. La parallélisation étant hors de portée d'un compilateur au- 
tomatique, elle doit être explicitement programmée, le "passage de message" 
s'avérant le mieux adapté. 
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